Exact Free Energy Functional for a Driven Diffusive Open Stationary Nonequilibrium 

System 



(N 
O 
O 
(N 



o 

a 

■i— > 
c3 



^3 

o 
o 



> 

m 
O 

<N 
O 
-i— > 



i 

C 

o 
o 



X 



B. DerridaQ 

Laboratoire de Physique Statistique, 24 rue Lhomond, 75231 Paris Cedex 05, France 

J. L. Lebowitzfj] and E. R. Speer| 

Department of Mathematics, Rutgers University, New Brunswick, NJ 08903 

We obtain the exact probability exp[— LJ r ({p(x)})] of finding a macroscopic density profile p(x) in 
the stationary nonequilibrium state of an open driven diffusive system, when the size of the system 
L — > oo. T , which plays the role of a nonequilibrium free energy, has a very different structure from 
that found in the purely diffusive case. As there, T is nonlocal, but the shocks and dynamic phase 
transitions of the driven system are reflected in non-convexity of T , in discontinuities in its second 
derivatives, and in non-Gaussian fluctuations in the steady state. 
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The behavior of macroscopic systems which carry 
steady currents is one of the central problems in nonequi- 
librium statistical mechanics Q. Of particular interest 
are stationary nonequilibrium states (SNS) maintained 
by contact with infinite reservoirs at the system bound- 
aries and subjected to bulk driving forces. A paradigm of 
such systems is a fluid in contact with a thermal reservoir 
at temperature T a at the top and XJ, > T a at the bottom, 
for which gravity supplies the bulk force (the Rayleigh- 
Benard system |Q); the system exhibits dynamic phase 
transitions corresponding to the formation of different 
patterns of heat and mass flow as the parameters are 
varied. By contrast, if T a > Tf, the system has no insta- 
bilities. These dynamic transitions are not understood, 
at present, in terms of a microscopically derived free en- 
ergy, despite various promising attempts ||. Here we 
obtain the analogue of such a free energy for the SNS of 
a model system which, despite its simplicity, has some 
dynamic transitions. 

We consider the one-dimensional asymmetric simple 
exclusion process (ASEP) on a lattice of L sites |). 
Each site i, i = 1, . . . , L, is either occupied by a single 
particle (r, = 1) or is empty (r-j = 0). In the interior of 
the system (2 < i < L — 1), a particle attempts to jump 
to its right neighboring site with rate 1 and to its left 
neighboring site with rate q (with < q < 1), succeeding 
if the target site is empty. At the boundary site i = 1 
(i = L) particles jump only to the right (left). These 
boundary sites are also connected to particle reservoirs: 
if site 1 is empty, it becomes occupied at rate a (by a 
particle from the left reservoir); similarly, if site L is oc- 
cupied, the particle may jump into the right reservoir at 
rate (3. 

This dynamics produces an SNS for which we calculate 
the large L behavior of Pl({p{x}) ~ exp[— LJ 7 ({p(x)})], 
the probability for observing microscopic configurations 
corresponding, in the limit L — > oo, x = i/L, to 
the macroscopic density profile p(x), < x < 1. 



J 7 ({p(x)}) is generally called the large deviation func- 
tional (LDF) in the mathematical literature ||; one al- 
ways has J 7 ({p(x)}) > 0, with equality holding only if 
p(x) — p(x), a typical density profile in the system, 
so that atypical profiles are observed with exponentially 
small probability for large L. For equilibrium systems 
T({p(x)}) is given as the difference in free energies for 
p{x) and for p(x). 

In our previous work we considered the case q = 1 , 
in which the bulk dynamics are symmetric. We obtained 
there an exact F s ({p}) which, unlike the free energy in 
equilibrium, was nonlocal, reflecting the generic presence 
of long range correlations in SNS || . Due to the purely 
diffusive nature of the bulk dynamics, J- s did not exhibit 
any phase transitions or instabilities. This is quite differ- 
ent from the asymmetric case considered here which has, 
due to the driven nature of the bulk dynamics, not only 
long range correlations but also a rich phase diagram in- 
cluding phase transitions and shocks ||, ||, [)). It is thus 
closer to a real fluid and gives rise to a correspondingly 
more complex T{{p\). 

Before describing our new results about T we sum- 
marize some known properties of the open ASEP [flo|| . 
We restrict ourselves in this paper to a — (1 — q)p a and 
(3 = (l — q)(l—pb) with < p a , pb < 1; the parameters p a 
and pt represent the densities in the left and right reser- 
voirs. For p a = pb = r the stationary measure is just 
a product measure [|[ [l^] with uniform density r. This 
means that all static (i.e. single time) properties of the 
system, including T , are the same as for an equilibrium 
lattice gas with fugacity z = r/(l — r). For this system 
the LDF is given by % 



Fea({p(x)}) 



p(x) loj 



p{x) 



+ {i-p{x))H 



1 - p(x) 



1-r 



(1) 
dx. 



Now recall that if an infinite system with ASEP dy- 
namics is started in an initial state with a macroscopic 
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where the sup is over all monotone nonincreasing func- 
tions F(x) which satisfy p a > F(x) > pb, < x < 1; for 

Pa < Pb, 



FIG. 1: The phase diagram of the open ASEP 



density profile p{x, 0) = p a for x < X, p(x,0) = pb for 
x > X, then when p a < pb the time evolved p(x,t) will 
maintain a sharp shock which will move with velocity 
V = (1 - q)(l - p a - Pb), while when p a > p b , p{x, t) will 
smooth out as in a "fan": p(x,t) = p a if x < x a (t), 

p{x,t) = p a + (p b - p a ){x - X a {t))/{x b {t) ~ X a (t)) if 

x a (t) < x < Xbit), and p{x,t) = pb if Xb{t) < x, with 
x a (t) = X + (1 — q)(l — 2p a )t, a = a,b. In each case 
(unless V = 0) the system will attain, as t — > oo, a 
constant density p in any hnite region. This gives an 
understanding || of the phase diagram shown in Fig. 1 
for the system with open boundaries. For p a < Pb, the 
shock will move to the right boundary when V > 0, leav- 
ing behind the constant density p — p a (phase Ai), and 
to the left boundary when V < 0, yielding p = pb (phase 
Bi). On the line S (V — 0) a typical p{x) is no longer 
constant, but corresponds to a shock at some point s, uni- 
formly distributed on [0, 1], where p(x) jumps from p a to 
Pb, i.e. p{x) = p s (x) = p a Q(s — x) + pbO(x — s) with 
the Heaviside function jnj. This line corresponds to a 
hrst order phase transition, with p discontinuous across 
S. For p a > Pb (phases A 2 , B 2 , and C) one sees the con- 
stant density, p a , Pb, or 1/2, which would have resulted 
from the fanlike behavior in the infinite system. 

We now turn to our new results. The line p a — pb, 
which separates what we will call the shock region p a < 
Pb from the fan region p a > Pb, is irrelevant for p but 
plays a crucial role in the LDF. Dchning 

g(h, f) = h\og[h(l - })] + (l-h) log[(l - h)f], (2) 
K(p a ,pb) =logp(l -p), (3) 

with p given in Fig. 1, we hnd: for p a > pb, 

T{{p[x)Ypa,Pb) (4) 

= -K{p a ,Pb)+ sup / dxg(p(x),F(x)), 

F(x) JO 



F{{p{x)};Pa,Pb)) = -K(p a ,Pb) 

+ W dxg(p{x),p a )+ / dxg(p{x),p b ) 

o<y<i [Jo J y 



(5) 



The fact that the function F(x) in (|]) is required to be 
monotone makes it, like y in dq), depend on the global 
form of p{x), so J 7 is a nonlocal functional of p(x). It 
can be shown |ll| that the optimal F in (|j), which we 
denote by F p , is constructed as follows: let G p be defined 
as the concave envelope of the function f Q (1 — p(y)) dy; 
then G' p (x) is monotone, 

F p {x) = G' p {x) if Pb < G' p (x) < p a , (6) 

and F p (x) = p a (F p (x) = p b ) where G' p (x) > p a 
(G' p (x) < pb). Note that F p (x) need not be continu- 
ous; it will generally consist of strictly decreasing pieces, 
where F p (x) = 1 — p(x), flat pieces, where F p (x) = p a , 
F p (x) — pb, or F p {x) is obtained from 1 — p(x) by 
a Maxwell construction (the integrals of F p {x) and of 
1 — p(x) over the latter intervals being equal), and pos- 
sible jumps downward. 

Before going on to discuss the derivation of ([|) and (||) 
we describe some consequences. 

(a) It can be verified that F({p{x)}\ p a , Pb) > 0; equal- 
ity occurs only when p(x) = p as given in the phase 
diagram, except on the first order line S, where it is 
the shock (typical) configurations p s (x) which satisfy 
T{{p s (x)}- Pa ,pb) = for all s S [0,1]. 

(b) J r ({p(x)}; p a , Pb) given by (^) is a convex functional 
of p(x) in the fan region p a > pb, since it is the supremum 
over F of J Q g(p,F)dx with g a locally convex function 
of p for every F. This is similar to what happens in 
the symmetric case 0. In the shock region p a < Pb, 
on the contrary, this is not true; for every p a ,Pb there 
are profiles p(x) near which T is not convex. This is 
easily verified on the line S where a superposition p(x) = 
Xp s (x) + (1 — X)pt(x), s ^ t, satisfies J 7 ({p(x)}) > for 
< A < 1. 

(c) The LDF in the fan region p a > Pb has similari- 
ties besides convexity to that in the symmetric case. In 
particular it is easy to see from (||) that T({p(x)}) > 
J r eq ({p(x)}) , where T eq is given by (1) with r replaced 
by the appropriate p. But in the shock region, p a < pb, 
this inequality is reversed: T({p(x)}) < T eq {{p{x)}) (as 
is clear from (||), since in region B\ (^4i), y = (y = 1) 
gives J-" cq ({p(a;)})). This is similar to what fluctuating 
hydrodynamics predicts for the Rayleigh-Bcnard prob- 
lem discussed earlier: fluctuations are decreased when 
T a > Tb (at least when T a — Tb is very small) and are 
increased when T a < Tb, even in the stable conductance 
regime p~2|. 
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(d) In the symmetric case discussed in |T^], as in an 
equilibrium system not at a phase transition (in any di- 
mension), the probability of small fluctuations about the 
typical state can be obtained from T as a limit. More 
precisely if we write p[x) — p(x) + -^=u(x) and then ex- 
pand T to second order (the first order term being zero) 
we get a Gaussian distribution for u(x) with covariancc 
C(x,x'), where C~ 1 (x,x') = 8 2 T '/ '8p(x)5p(x') evaluated 
at p = p. This covariance is the suitably scaled micro- 
scopic truncated pair correlation B . For the asymmetric 
case discussed here the distribution of fluctuations need 
no longer be given by the LDF; in fact, 8 2 J-/8p(x)8p(x') 
is discontinuous at p = 1/2 in the interior of region C of 
the phase diagram, i.e., where p a > 1/2 > pb- Further- 
more, the fluctuations there are no longer Gaussian. 

To see this discontinuity in an explicit example let 
p(x) — 7} + eQ(x — s), with pb < \ — e < p a ; here it is easy 
to compute T(e) = J r ({p(x)}). First, F p {x) = 1 — p(x) if 
e > and F p (x) — k — e(l — s) if e < 0; the constancy of 
Fp for e < is due to the concave envelope construction 
(0). Then from (0) (we give only the small-e behavior): 



P{e) 



4(1 - s)e 2 + . 
4(1 - a) (1 - 



if e > 0, 
if e < 0. 



(7) 



The discontinuity of d 2 !F(e) / de 2 at e = signals that the 
fluctuations are anomalous (non-Gaussian). (Note that 
for s = 0, d 2 !F(e)/de 2 is continuous at e = and is in fact 
equal to the inverse of the variance in the total number 
of particles |pTf.) 

These non-Gaussian fluctuations can be observed by 
considering the total number N s of particles on lattice 
sites in [s, 1], i.e., sites i with sL < i < L. A calculation 
using the results of JTJ| for the microscopic probabilities 
shows that [N s - (1 - s)£/2]/\/I converges, as L — > oo, 
to a random variable p with a well defined but nonsym- 
metric (and non-Gaussian) distribution having density 



s3/2(x- a )s 



dtt 2 e 



(8) 



This is in contrast with equilibrium systems (not at a 
phase transition) for which statistical mechanics |L4j pre- 
dicts Gaussian fluctuations of the number of particles in 
a macroscopic region (the variance being related to the 
compressibility). For large values of \p\, (||) yields 



\ogp(p) 



4p 2 
1-s 



Ap 2 



(1 -*)(! + *) 



if M> 1, 
if/i<-l. 



(9) 



This agrees with the results of a large deviation calcu- 
lation |ll[ of the probability exp[— LT(e)] of observing 
mean density ^ + e in the interval [s, 1], with no other 
constraints: the small-e behavior of jF(e) agrees with 
the large-/Lt asymptotics of (^|) under the identification 



p = \/L{l — s)e. (Note that (||) differs from (Q) because 
for (Q) a constraint is imposed also in the region [0, s\.) 

Derivation: To obtain T in (|]) and (^) we use the 
exact expression for the measure Pl{t) provided by the 
matrix method However, rather than calculating the 
probability of a given macroscopic profile p{x) directly by 
summing Pl(t) over all configurations r corresponding 
to that profile, as we did for the symmetric case in we 
follow here a different path, which has its origin in an a 
posteriori observation made in |Q. We noted there that 
while J 7 s ({p(x)}; p a , Pb) is nonlocal, it possesses a certain 
"additivity" property which, if it could have been es- 
tablished independently, would have yielded T s . This is 
exactly what we do for the ASEP: we first derive an addi- 
tivity property from the matrix representation of Pl (r) , 
then use the additivity to obtain T . Full details are given 
in Hj; here we give only a partial sketch of the argu- 
ments. 

Let us introduce the LDF T\ a w ({p(x)}; p a , pb) — 
L" 1 log P L ( b -a)({p{x)}; Pa, Pb) for a system of L(b - a) 
lattice sites in contact with reservoirs at densities p a and 
Pb- Let K(p ai pb) be as in (||) and define 

ft[a,b]({p(x)}; Pa, Pb) (10) 
= F [a , b ]({p(x)}; p a , Pb) + (b- a)K(p a ,Pb)- 

It is shown in Jll| that for p a > pb we have 



T-(-[a,b]{{p(x)};Pa,Pb) 



sup 

Pb<Pa<Pa 



n[a,c]({p(x)}; Pa, Pc) 
11) 



+ K[c,b]({p(x)}; Pc, Pb 
while for p a < p b , 

H[a,b]({p{x)};Pa,Pb) (12) 

= min [U [a . c] {{p{x)}; p a , Pa) + H [c M]({p(x)}; p a , Pb), 

T-(-[a,c] {{p(x)Y, Pa, Pb) + H[ c ,b] {{p(x)Y, Pb, Pb)] ■ 

Equations ([T^) and ( p^ ) are the additivity relations for 
the ASEP. 



In this letter we will sketch a derivation of (11) in the 



special case q = 0; for the general case and for the deriva- 
tion of (12) see |Q. We use the matrix formula || for a 
system of N sites: 



Pn(t) 



(Painl^Dn + Eii-r^lpb) 



(p a \(D + E)K\ Pb } 
where the operators D, E and vectors {p\, \p) satisfy 



DE = D + E, 



(p\E=-( P \, 
P 



D\p) 



1-P 



\P)- 



(13) 

(14) 
(15) 



If ( J15| ) is extended to complex values of p we may write 
the exact additivity formula 



(palXoX^Pb) 
(Pa\Pb) 



(16) 
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1 

2tH 



{Pa ~ Pb)dp (p a \X \p) (pjX^Pb) 
(Pa-P)(p-Pb) (Pa\p) (p\Pb) 



where Xq, X\ are arbitrary polynomials in D and E and 
the contour is a circle \p\ = R with p\, < R < p a . To 
obtain ([l6|), note that it suffices to take X t — E Pi D Qi , 
since from ( |i~4"| ) any polynomial in D and E can be written 
as a sum of such terms. The cases go = or p\ = are 
immediately obtained from ( p"5| ) and the residue calculus; 
the general case follows by an inductive argument, as the 
case qo,Pi can be reduced to the cases qo — l,pi and 
go, Pi — 1 using ( (Til) on the left hand side of (|l^ ) and the 
corresponding identity (1 — p)~ 1 p~ 1 = (1 — / o)~ 1 +p~ 1 on 
the right hand side. Now the weights (p a \X \p) / (p a \p) 
and (p\Xi\pb)/ (p\pb) in (|l6| ) are polynomials in 1/(1 — p) 
and 1/p, respectively, with positive coefficients, so that 
the integrand has a Laurent series for pb < \p\ < p a 
with positive coefficients and hence is a convex function 
of real p for Pb < p < p a \ the minimum p mm for such real 
p (which must occur since there are poles at p = p a and 
p = pb) is a saddle point for the complex integral. Thus, if 
p(x) is a given profile and Xq (respectively X\) is a sum 
of products of L(c — a) (respectively L(b — c)) factors 
of D or E consistent with the left (respectively right) 
part of this profile, and we assume that that the weights 
depend exponentially on L, we expect the integral to be 
dominated by this saddle point, leading to 



i_ lo {p a \xMpb) 

L (Pa\pb) 



(17) 



inf 

Pb<P<Pa 



1 (Pa\X Q \p) 
L (Pa\p) 



l lQg 

L {p\pb) 



Using (|l7|), (|l3|), and the relation 
( Po |(I> + ^ 



TV" 1 log ■ 



(Pa\Pb) 



-K(p ai p b ), 



(18) 



we obtain (|ll]). Thus the rather surprising supremum in 
( |ll| ) , which corresponds to a choice of the least probable 
alternative, arises mathematically through the contour 
representation (16) and its domination by a real saddle 
point; we still lack a physical explanation of (|l7j). 

To go from (|ll|) to (j|) we divide our system into n parts 
of equal length, the k th interval being [(k — l)/n,k/n], 
k = 1, . . . , n. Now note that for very large n most of these 
intervals must have reservoir densities pk-i,Pk which are 
nearly equal, and that the LDF for these intervals is ap- 
proximately given by (|l|) with r ~ Pk-i — Pk- On the 
other hand, the total length of the intervals for which 
this is not true (corresponding to points of discontinuity 
in the function F p ) will approach as n — > oo; these con- 
siderations lead then directly to (Q) in this limit. We pass 
from (12) to (|5|) by a similar process of subdivision, but 
the argument is even simpler since each interval except 
one has equal reservoir densities. 

Conclusion: Tt would he interesting to know how key 



for the LD functions, the suppression or enhancement 
of deviations of the density profile from its typical value 
as the reservoirs and the internal field cooperate or com- 
pete, and the non-Gaussian fluctuation (||) of the number 
of particles in a box of length cL, < c < 1 — extend to 
more realistic systems, and whether they could be un- 
derstood by a dynamical approach, as is done for the 
symmetric case in Jll| 
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